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Abstract 

Factorization  of  linear  programming  (LP)  models  enables  a  large  portion  of  the  LP 
tableau  to  be  represented  implicitly  and  generated  from  the  remaining  explicit  pari.  Dy¬ 
namic  factorization  admits  algebraic  elements  which  change  in  dimension  during  the  course 
of  solution.  A  unifying  mathematical  framework  for  dynamic  row  factorization  is  presented 
with  three  algorithms  which  derive  from  different  LP  model  row  structures:  generalized 
upper  bound  rows,  pure  network  rows,  and  generalized  network  rows.  Each  of  these  struc¬ 
tures  is  a  generalization  of  its  predecessors,  and  each  corresponding  algorithm  exhibits  just 
enough  additional  richness  to  accommodate  the  structure  at  hand  within  the  unified  frame¬ 
work.  Implementation  and  computational  results  are  presented  for  a  variety  of  real-world 
models.  These  results  suggest  that  each  of  these  algorithms  is  superior  to  the  traditional, 
non- factorized  approach,  with  the  degree  of  improvement  depending  upon  the  size  and  qual¬ 
ity  of  the  row  factorization  identified. 
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1  Introduction 


A  recurring  theme  in  the  development  of  algorithms  for  linear  programming  has  been  the  iden¬ 
tification  and  exploitation  of  special  problem  structure,  ideas  as  apparently  disparate  as  the 
bounded- variable  simplex  method,  primal  and  dual  decornpo-sition  methods,  pure  and  gener¬ 
alized  network  primal  simplex  algorithms,  primal  partitioning  and  column  generation  schem(‘s 
may  be  unified  to  a  degree  with  this  view. 

The  factorization  approach  introduced  by  Craves  and  McBride  [1976]  Isolates  sjrecial  .struc¬ 
ture  in  LP  tableaus.  We  are  interested  in  using  factorization  to  reinterpret  existing  algorithms, 
and  to  discover  common  principles  and  apply  them  to  develop  new  algorithms.  Although  all 
algorithms  developed  tais  way  will,  in  theory,  solve  any  LF^,  the  efficiency  of  any  particular  fac¬ 
torization  approach  will  be  influenced  by  the  relative  number  of  factored  constraints  and  their 
influence  on  the  algorithm;  the  .‘lize  and  quality  of  the  special  structure  i.solated  determines 
the  influence  of  any  particular  factorization  applied  to  any  particular  LF’. 

Based  on  prior  work  by  Brown  and  Craves  [1975],  in  which  generalized  upper  Iround  rows 
were  successfully  incorporated  in  a  large-scale  optimization  system,  we  are  interrs^ted  in  pursuing 
dynamic  row  factorization,  where  the  dimension  of  the  factored  structure  may  vary  (or  even 
fail  to  be  present)  as  the  solution  progresses.  In  our  setting,  we  require  the  row  structure  of 
the  model  instance  to  be  specified  prior  to  .solution,  and  that  this  structure  remain  fixed  during 
solution.  An  extension  of  this  approach  is  to  allow  the  row  structure  to  vary  as  the  model  is 
solved;  this  is  a  conceptually  simple  extension  of  the  approach. 

Each  algorithm  is  developed  by  factoring  the  constraints  of  the  LF’  model  into  two  classes; 
those  that  have  the  special  structure  (factored)  and  those  that  do  not  (eTplic.it).  This  constraint 
factorization  induces  a  factored  structure  in  the  LP  tableaus  which  is  exploited  computationally. 
We  demonstrate  the  dynamic  factorization  approach  for  three  special  structures: 

generalized  upper  bound  rows; 

pure  network  rows;  and 

generalized  network  rows. 

We  implement  each  of  the  factorization  algorithms  by  integrating  it  within  the  X-Systein 
(Brown  and  Craves[1975]). 

While  the  terms  “partitioning”  and  “factorization”  are  frequently  used  interchangeably  in 
the  literature,  we  observe  a  distinction  between  the  two  approaches.  We  consider  partitioning 
methods  to  be  based  on  special  structure  in  the  original  problem  instance,  which  need  not 
induce  special  structure  in  the  LP  tableau — in  fact,  the  method  need  not  be  taljleau- based. 
In  contrast,  factorization  methods  are  based  on  special  structure  which  occur";  in  bases  and 
thus  in  the  basic  tableau.  Thus,  we  clas.sify  dual  decomposition  (Dantzig  rad  Wolfe  [1960]), 
primal  decomposition  (Benders  [1962]),  and  primal  partitioning  (Rosen  [1964])  as  examples  of 
partitioning  methods. 

Perhaps  the  earliest  example  of  what  we  consider  factorization  is  the  treatment  of  simple 
upper  bounds  by  Dantzig  [1954]  and  [1963]  and,  independently,  by  Charnes  and  Lcmke  [1954]. 
They  observe  that  it  is  more  efficient  to  enforce  the  “logical”  upper  bound  constraints  with 
logical  tests  within  the  algorithm  rather  than  treat  them  explicitly  along  with  other  “structural” 
constraints.  While  not  originally  presented  in  the  context  of  a  formal  tableau  factorization,  the 
approach  is  easily  viewed  as  such. 

The  mutual  primal-dual  method  of  Clraves  [1965]  focuses  attention  on  t  he  special  role  of 
nonnegativity  constraints  in  linear  programming.  A  clear  distinction  is  drawn  betw'een  the 
computational  convenience  of  treating  nonnegativity  constraints  implicitly  rather  than  explic¬ 
itly  and  the  unambiguous  mathematical  equivalence  of  all  prol>lem  constraints,  structural  or 
nonnegativity.  Emphasizing  the  special  importance  of  inetjuality  constraints,  the  approach 
yields  an  elegant  theory  and,  as  we  will  .see,  efficient  implementations.  We  view  this  algorit  hm 
as  the  first  formal  example  of  factorization. 
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A  similar  primal-dual  algorithm  is  presented  by  Balinski  and  (Jomory  [1965].  Relatwi  work, 
in  which  efforts  are  made  to  exclude  slacks  from  the  product-form  representation  of  the  primal 
basis,  includes  that  of  Zoutendijk  [1970]  and  Powell  [1975]. 

Dantzig  and  Van  Slyke  [1967]  extend  the  earlier  work  for  simple  upper  bounds  and  lend 
a  more  structured  treatment  to  generalized  upper  bounds  (Cll'B).  In  a  problem  with  p  CUll 
constraints  and  m  structural  constraints,  their  approach  requires  a  working  liasis  of  dimen.sion 
(m  -f  1),  a  considerable  savings  when  p  is  large. 

Hartman  and  Lasdon  [1972]  specialize  this  approach  to  the  multicommodity  capacitate*! 
transshipment  problem.  In  this  case,  the  structure  of  the  basic  pure  network  columns  in¬ 
troduces  additional  structure  into  the  working  basis,  allowing  further  simplifications  in  basL< 
representation  and  update  techniques.  Helgason  and  Kennington  [1977]  develop  technicjues  for 
representing  the  working  basis  in  product  form  and  provide  graphic  interpretation  of  the  basis 
updates.  Kennington  [1977]  reports  an  implementation  of  the  algorithm. 

McBride  [1972]  and  Graves  and  McBride  [1976]  formalize  and  generalize  the  factorization 
approach.  They  view  factorization  as  a  unifying  framework  for  tableau- based  simplex  spe<  ial- 
izations  and  illustrate  this  by  developing  a  variation  of  the  GUB  algorithm  of  Dantzig  and  Van 
Slyke  and  a  GUB  algorithm  for  the  doubly-coupled  linear  programs  of  Hartman  and  Lasdon 
[1970].  They  present  a  new  algorithm  for  the  .set  partitioning  LP  and  an  equality-constrainrsi 
form  of  the  pure  network  with  side  constraints  model.  Brown  zind  Graves  [1975]  report  an 
implementation  of  inequzdity-form,  dynamic  GUB  row  factorization  for  large-scale  problems. 

Schrage  [1975]  extends  the  succession  of  simple  and  generalized  upper  bounds  by  introducing 
varirble  upper  bounds  (VUB),  which  are  coiLstraints  of  the  form  <  x*,  where  x*  is  said  to  be 
the  variable  upper  bound  of  Xj.  Schrage  implicitly  represents  the  VUB  constraints  by  expressing 
VUB  variables  in  terms  of  other  variables.  This  permits  the  basis  representation  to  be  treated  in 
two  parts,  one  a  large  matrix  which  changes  infrequently  and  thus  needs  only  occasional  update, 
and  the  other  a  small  working  basis  which  requires  regular  attention.  Thus,  computation  and 
storage  savings  may  be  realized.  Schrage  [1978]  extends  these  ideas  to  what  he  calls  generalized 
VUB  (GVUB)  constraints,  which  arise  frequently  in  models  with  fixed  charges. 

Klingman  and  Ru.s.sell  [1975]  sketch  a  factorization  method  for  solving  transportation  prol>- 
lems  with  side  constrmnts.  They  suggest  techniques  for  performing  simplex  iterations  and 
updating  the  problem  representation.  Chen  and  Saigal  [1977]  present  a  similar  approach  for 
solving  capacitated  network  flow  problems  with  additional  linear  constraints.  Both  of  the.se  pre- 
.sentations  employ  a  graphical  description  of  the  basis  update  and  treat  the  basis  in  two  parts: 
one  corresponding  to  a  rooted  spanning  tree  defined  on  the  underlying  graph,  and  the  other 
a  general  working  basis.  Glover,  et  al.  [1978]  report  an  implementation  of  the  Klingman  and 
Rmssell  design,  but  one  which  (curiou-sly)  only  accommodates  a  single  side  constraint.  McBride 
[1989]  reports  an  implementation  which  requires  the  pure  network  rows  to  be  equalities  and 
allows  more  than  one  side  constraint. 

Generalized  networks  with  .side  constraints  are  addres.sed  by  Hultz  and  Klingman  [1976], 
who  present  details  for  the  simplex  priceout,  column  generation,  and  basis  update.  Hultz  and 
Klingman  [1978]  report  an  implementation  that  (curiously)  .solves  the  “singularly  constrained” 
generalized  network  problem.  McBride  [1989]  reports  an  implementation  that  is  not  restricted 
to  a  single  side  constraint. 

The  factorization  approach  has  been  extended  by  consideration  of  embedded  structures. 
Glover  and  Klingman  [1981]  consider  an  LP  with  embedded  pure  network  structure,  i.e.,  the 
pure  network  structure  appears  in  only  a  subset  of  the  rows  and  columns  of  the  technological 
coefficient  matrix.  They  give  an  algorithm  similar  in  spirit  to  their  pure  netw'ork  with  side 
constraint  model,  but  the  presence  of  the  “side  variables”  significantly  complicates  the  basis 
representation  and  update.  They  report  an  implementation  of  the  algorithm  but  (curiously) 
restrict  test  problems  to  have  no  complicating  variables. 

McBride  [1985]  solves  an  LF’  with  enibedded  generalized  network  structure,  presenting  meth¬ 
ods  for  pricing,  column  generation,  basis  representation  update  and  data  struct  ures.  A  .succes.s- 
ful  implementation  is  reported  to  be  about  five  times  faster  than  MINO.S  (ca.  1977:  Murtagh 
and  .Saunders  [1977])  for  the  models  tested. 
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Algorithms  to  solve  problems  with  special  sulistructures  have  molivatwl  r<^earch  to  effi¬ 
ciently  identify  such  substructures.  Brearley,  Mitra  and  Williams  [I975j  descril)e  algorithms 
for  detecting  CJUB  row  sets  and  exclusive^row’  structure  sets  (a  set  of  rows  whose  structure 
may  be  transformed  to  GUB  by  column  scaling).  Greenberg  and  Rarick  []974j  and  Brown  and 
Thomen  [1980]  develop  algorithms  to  identify  GUB  sets.  Brown  and  Wright  [198-1]  identify 
pure  network  constraint  substructures.  Brown,  McBride  and  Wood  [1985]  present  a  method  for 
locating  embedded  and  row-only  generalized  network  structures. 

Todd  [1983]  develops  a  geometric  interpretation  of  factorization  which  Ls  for  our  purpo.ses 
equivalent  to  the  algebraic  development  of  Graves  and  McBride  [1976]. 

In  the  following  .sections  we  establish  notational  conventions,  develop  the  mathematical 
foundation  of  a  primal-dual  simplex  method  and  show  the  effects  of  row- factorization.  Next, 
a  general  row-factorization  algorithm  is  developed,  specialized  to  GliR,  pure  network,  and 
generalized  network  rows,  and  tested  on  a  suite  of  real-world  problems.  Finally,  we  discuss  how 
the  methods  can  be  generalized  further  and  how  they  can  be  applied  more  eff»?ctively. 


2  Mathematical  Preliminaries 

The  traditional  statement  of  the  linear  programming  (LP)  prol)lem  is 


(Lf’)  min  :  wy 
y 

s.t.  Oiy  <  r'l  .  f  =  l . m 

Cjy  >  0  .  j  =  1 . Tj  . 

where  y  is  an  n-vector  of  decision  variables,  w  a  vector  of  cost  coefficients,  each  a,  an  n-vec  tor 
of  technological  transformation  coefficients,  each  rj  a  scalar  right-hand  side  coefficient,  and  r.j 
the  unit  vector.  While  this  statement  of  the  problem  is  clear  and  unambiguous,  there  arc 
reasons  for  preferring  an  alternative.  The  insistence  upon  drawing  a  formal  distinction  between 
the  “structural”  constraints  Ojy  <  Tj  and  the  “non negativity”  constraints  e.^y  >  0  oljscures 
the  mathematical  structure  of  the  problem  by  suggesting  that  the  two  types  of  constraints  are 
inherently  different.  Certainly  the  exploitation  of  the  special  structure  of  the  e^y  >  0  constraints 
leads  to  computational  efficiencies;  however,  in  our  theoretical  development  of  the  algorithm, 
we  prefer  to  treat  them  simply  as  general  inequality  constraints. 

In  order  to  achieve  a  consistent  form,  we  rewrite  the  nonnegativity  constraints  as  -  ^■jy  <  0 
and  group  them  with  the  structural  constraints.  The  problem  statement  then  becomes 


(PLP)  min  ;  wy 

s.t.  :  a,y  <  n 


,  i  =  1 . m  +  n. 


where  wy  is  called  the  extremal  function. 

From  the  standpoint  of  a  primal  algorithm,  a  matrix  partitioned  form  of  the  primal  tableau 
is  derived.  Let  {a,|,aij, . . .  ,a,^}  be  a  basis  for  R”  at  j/®.  For  notational  convenience  we  will 
partition  the  constraints  into  two  .sets,  those  that  are  basic  (binding)  at  y”  and  those  that  are 
nonbasic  (not  necessarily  binding)  at  y® 


■ 
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D  = 

•«  ■« 

;  9  = 
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zz: 

'  u„., 

Qm 

Using  this  notation,  the  current  basic  solution  i/°  may  be  expressed  as  Bip  =  /,  and  since 
the  rows  of  B  are  by  definition  linearly  independent,  j/°  exists. 

To  isolate  the  important  algebraic  components  let  us  assume  that  at  the  current  i)asir 
solution  the  basis  consists  of  h  structural  constraints  and  (n  —  h)  nonnegativity  constraints. 
Then 
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bh 
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In  partitioned  matrix  form, 


h 


n~h 


Ai2  I 

0  -I  )  }n~h  . 


and  thus 

}h 

}n  -  h 

Similarly,  D  can  be  written  as 


p  =  = 


n~h 


V  0  / 
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51 
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and  in  partitioned  matrix  form 


Then 


D  = 


h  n~h 


( 
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I  0 

V  -^21  ^7  ) 


]h 

}  rn  -  h  . 


DP 


■DB'  = 


4-' 

^11 


a;}  A 


n  ^12 
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A22  -  A2iA[i’A\2  J  }  ni  —  h 


and  we  shall  call  DP  the  principal  part  of  the  tableau.  By  partitioning  in  =  (u’j.uii).  — 
and  =  (y?,y2)'  ^he  complete  tableau  may  be  written  in  partitioned  matrix  form  as 


h 

( 

n  —  h 

1 

\ 

A11A12 

y? 

-A2\AjI 

A22  -  A2iAJ','Ai2 

52  -  A2iy?  -  ^2252 

}  m  —  /? 

V 

W2  -  W]Ai'Ai2 

1 

0 

}1  • 

Note  that  y®  is  displayed  explicitly  in  this  tableau.  Also,  y^  =  0  since  the  corresponding 
nonnegativity  constraints  are  basic  and  thus  binding. 

The  corresponding  dual  problem  is 


(DLP)  max  :  xr 

s.t.  xa^  <  ,  j  =  ,n 

xe,  <0  ,  i  =  l,...,m. 

To  develop  a  matrix  partitioned  form  of  the  dual  tableau,  we  proceed  as  before.  Assuming 
the  dual  basis  consists  of  h  structural  constraints  and  m  —  h  nonnegativity  constraints,  we  have 


T  . /m)  =  (a2•,a-'^...,  a",  . e,'") 

u  =  (u*.  . . . ,  . 0), 


u  =  (u\  u^)  =  (ui',  0). 

The  nonbasic  constraints  are  then 

K  =  . fc")  =  (c",e‘’ _ _ n^") 

V  =  (vKv^ . u")  (0,0, - 0,  lil-"'*!, _ ll,>2’'). 
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so  I!  =  ,  v^)  =  (0,  w^)  . 

The  matrix  partitioned  form  of  the  basis  is  then 


T  = 


/in 

/iai 


0 

f  / 


}h  . 

}  TTl  —  h  , 


and  with  the  choice  of  Q  =  T  ' 


h  m-h 


\-A2xA;,'  I  I 


with  the  remaining  constraints  forming 


}m  —  h  . 


[  A 
0  A22 


The  principal  part  of  the  dual  tableau  is 


QK 


! 

0 

_ 1 

■/  >1,2  ■ 

-A2\Aii  I 

0  i422 

/in'  /i]l'/il2 


[  ~i42l/ii|'  /i22  -  /i2l/ill'/il2  J  ' 

which  we  find  to  be  exactly  the  principal  part  of  the  primal  tableau,  so 

DP  =  QK  . 

Thus  a  single  tableau  representation  supports  both  primal  and  dual  algorithms. 


3  Interpretation  of  Primal  and  Dual  Forms 

We  may  interpret  a  primal  or  dual  algorithm  as  simply  different  perspectives  of  this  same 
tableau,  wherein  a  primal  algorithm  basis  change  is  viewed  as  exchanging  primal  constraints  and 
a  dual  basis  change  exchanges  dual  constraints.  The  classical  Simplex  Method  may  then 
be  interpreted  as  solving  (PLP)  using  the  dual  perspective.  That  the  classical  Simplex 
Method  is  naturally  interpreted  as  a  dual  algorithm  comes  as  a  surprise  to  the  conventionally 
trained.  However,  the  consequent  mathematical  insight  is  compelling,  especially  in  light  of 
the  notational  simplification  and  apparent  underlying  role  of  which  we  refer  to  as  the 
transjorviation  kernel. 

There  are  several  reasons  for  preferring  a  Primal-Dual  algorithm  to  the  Simplex  Method. 
From  a  computational  standpoint,  because  slack  variables  are  carried  logically  rather  than 
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introduced  explicitly,  we  are  able  to  clearly  identify  the  essential  information  nef?ded  lo  exe<  ute 
the  algorithm.  The  matrix  plays  a  key  role  in  the  calculation  of  the  tableau,  and  the  entire 
tableau  can  be  constructed  from  i4,,*  and  original  problem  data.  Since  ylj,®  is  a  submalrix 
of  the  inverse,  T“*,  used  by  the  Simplex  Method,  it  is  smaller  and  recjuires  fewer  arithmetic 
operations  to  update  than  does  T~*. 

A  second  advantage  of  a  Primal-Dual  Algorithm  lies  in  the  flexibility  it  offers  for  special¬ 
ization  to  particular  problem  classes  or  structures.  Indeed,  it  is  the  special  structure  and 
simplicity  of  the  nonnegativity  constraints  that  motivate  the  development  of  the  algorithm  in 
the  first  place.  It  is  frequently  the  case  that  other  special  structures  can  be  identified  in  classes 
of  (PLP).  Examples  of  such  structures  include  simple  upper  bounds,  generalized  upper  bounds, 
variable  upper  bounds,  pure  and  generalized  network  substructures,  etc,  .Such  struc  ture  may  be 
“static”  in  that  its  nature  and  dimension  remains  fixed  throughout  the  solution  process,  or  the 
structure  may  be  “dynamic”  in  which  case  its  precise  nature  and/or  dimension  may  vary  as  the 
problem  is  solved.  Some  special  structures  may  be  more  strongly  characterized  by  their  c  (dumn 
structure  and  others  by  their  row  structure.  The  Primal-Dual  perspective  leads  naturally  to 
explanations  of  the  implications  of  virtually  any  such  problem  structure  and  greatly  simplifies 
the  implementation  of  such  a  specialization. 

When  an  LP  appears  as  a  subproblem  in  a  more  sophisticated  solution  setting  (for  examide, 
in  a  mixed  integer  programming  problem  or  a  nonlinear  programming  problem),  the  row/column 
symmetry  of  a  Primal-Dual  Algorithm  is  of  critical  importance  in  specializing  the  solution 
approach.  The  inherent  .symmetry  of  such  an  algorithm  permits  easy  adaptation  to  branch- 
and-bound  and  cutting-plane  approaches  to  mixed  integer  programming,  to  column  generation 
settings,  as  well  as  to  primal  and  dual  decomposition  techniques. 

We  believe  the  reason  for  this  flexibility  offered  by  the  algorithm  lies  in  its  more  complete 
mathematical  foundation.  There  is  a  natural  consistency  that  arises  from  the  choice  of 
a  vector  space  having  the  same  dimension  as  the  problem  variables  that  is  lacking 
in  other  approaches.  A  natural  geometric  interpretation  of  the  .solution  trajectory  follows  di¬ 
rectly  from  this  development.  Incidental  issues  such  as  finding  an  initial  basic  feasible 
solution  and  dealing  with  degeneracy  are  resolved  constructively  in  this  math¬ 
ematical  framework  [Graves  1965).  Other  approaches  re.x>rt  to  unnecessarily  complicate<i 
tangential  efforts. 

All  the  research  results  reported  here  can  be  developed,  with  some  effort,  in  the  framework 
of  the  classical  Simplex  Method.  However,  we  choose  to  pre.sent  these  results  in  the  manner  of 
their  development  —  the  mutual  Primal-Dual  view  presented  by  Graves. 


4  Column  and  Row  Generation 

Rather  than  maintain  a  complete  tableau  DP,  now  consider  the  generation  of  Just  column  c 
of  this  tableau.  Rewriting  in  a  manner  that  highlights  our  intentions,  and  labelling  row  and 
column  partitions  for  identification 

0)  Uj) 

(  \ 

[)P-  (0  [^n'l 

(if)  \ -A2i[A,,']  A22  -  A2i[A||' A12);  ■ 

By  properly  sequencing  our  computations  we  will  exploit  the  fact  that  region  (ii)  of  a  given 
column  is  simply  a  linear  combination  of  terms  in  region  (i)  of  the  same  column. 

Assume  we  want  to  place  the  current  representation  of  column  r  into  a  work  array  r,  which 
we  partition  as  z’^  =  (z,  to  correspond  to  regions  (i)  and  {ii).  Expre.s.sed  in  terms  of  the 
trEinsformation  kernel  A|“j',  we  compute  column  c  as 
if  c  is  in  (j), 


H 


and  then 


=  Mn’ r . 


-2  =  -A2i[::i]  ; 


or,  if  0  is  in  (jj), 


and  then 


-1  =  (^u  (^1:2)'' 


z-i  —  {M2Y  -  • 

Then  the  current  representation  of  column  c  is  available  in  =  (rf ,  zj). 

The  computation  of  7'ow  r  of  the  tableau  proceeds  in  a  similar  manner.  We  now  view  the 
principal  part  of  the  tableau  as 


DP  = 


_  (0 


0) 

0}) 

[Ail] 

[^h’M 

12 


(ii)  \[— i421i4,]]  2^22  +  [—^21^11  1"^!2  / 


If  we  want  to  place  the  current  representation  of  row  r  in  a  work  array  partitioned 
conformably  with  O')  and  (jj)  as  z  —  (23,  ^4),  we  compute 
if  row  r  is  in  (i), 


and  then 


or,  if  row  r  is  in  {ii), 


and 


^4  =  (•^3]>ii2  ; 
%  =  [(-v42l)ri4^,‘] 


Zi  =  {A22)r  +  [231^12  , 

and  the  current  representation  of  row  r  is  available  in  i  =  (£3,  £4). 

We  see  that  in  each  case  calculations  proceed  by  first  using  a  representation  of  /In'  to 
compute  a  portion  of  the  row  or  column  and  then  using  this  initial  computation  and  original 
problem  data  to  compute  the  remaining  part.  We  will  discover  that  our  specializations  ex¬ 
tend  this  approach  by  introducing  additional  tableau  partitions  which  allow  this  computational 
strategy  to  be  applied  on  a  larger  .scale. 
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5  Transformation  Kernel  Update 

The  dynamic  behavior  of  /t  j  i’  is  important.  VVe  see  from  tiie  primal  row  l)asis  H  an<i  nonbasif 
rows  D  that  the  dimension  of  v4,|'  corresponds  to  the  number  of  liasic  slnictura!  constraints, 
or,  equivalently,  to  the  number  of  nonbasic  nonnegativity  constraints  (recall  that  if  a  nonnega¬ 
tivity  constraint  is  nonbasic  and  thus  nonbinding,  the  corresponding  variable  may  i)o,ssibly  be 
nonzero).  Recalling  that  our  primal  view  of  a  basis  exchange  is  as  an  exchangf*  of  < onstraint.s 
between  B  and  D,  we  see  that  one  of  four  cases  may  (tccur  during  a  pivot 

A  structural  constraint  enters  the  basis  B  and  a  structural  constraint  leave's  the  basis  and 
enters  D.  Since  the  number  of  basic  structural  constraints  (and  the  numlter  of  /njnbasic 
nonnegativily  constraints)  remains  unchanged,  the  dimension  of  A,]'  is  unchanged.  A 
pivot  of  this  type  involves  a  row  in  region  (j)  of  B  and  a  row  in  regittn  (b)  of  I),  and  thus 
it  corresponds  to  a  pivot  coordinate  in  the  location  ((ti).(j))  of  the  tableau  DP. 

A  nonnegativity  constraint  enters  the  basis  and  a  nonnegativify  constraint  leaves  the 
basis.  Again,  the  dimension  of  Aj"!'  remains  unchanged.  Since  this  pivot  ttnolves  a  row 
in  region  {jj)  of  B  and  a  row  in  region  (t)  of  D,  the  corresponding  tableau  DP  pivot 
coordinate  lies  in  ((i),  (jj))- 

A  ,structural  constraint  enters  the  basis  and  a  nonnegativity  constraint  leavf*s  the  basis, 
and  thus  the  number  of  basic  .structural  constraints  (et^uivalenf  ly.  the  number  of  nonba'^ic 
nonnegativity  constraints)  uici'eases  by  one.  The  dimension  of  A|  |'  is  increased  by  one. 
This  corresponds  to  a  pivot  coordinate  in  region  ((ti).  {jj))  the  tal)!eau  DP. 

A  nonnegativity  constraint  enters  the  basis  and  a  structural  constraint  leavrs  the  basis, 
and  thus  the  dimension  of  Aj”,'  decreases  by  one.  The  corresponding  pivot  coordinate  in 

DP\s{{i)dj))- 

We  see  that  we  may  exert  some  influence  on  the  behavior  of  the  dimension  of  A,/  by  our 
strategy  for  selecting  target  exchanges  for  primal  and  dual  coastraints  (i.e..  our  oncing  strategy) 
and  through  our  tie-breaking  rules  for  choosing  pivot  row/column,  and  that  this  dynamism  is 
an  inherent  feature  of  an  effective  algorithm.  We  have  already  seen  the  fundamental  importance 
of  the  kernel  (An)  in  our  computations.  Thus,  a  .successful  implementation  must  manage  this 
dynamic  behavior  efficiently  and  reliably. 


6  Factorization 

The  row-factorized  problem  to  be  considered  is 


(FLP)  mm  : 

wy 

s.t.  : 

Ey  <  r 

}  explicit  constraints 

Ey  <  b 

}  factored  const  raints 

-ly  <  0 

)  nonnegativity  constraints 

where  y  is  an  n-vector  of  decision  variables,  w  a  vector  of  cost  coefficients,  E  a  matrix  of  con¬ 
straint  coefficients  for  “explicit”  constraints  with  right-hand  side  m-vector  r,  F  a  matrix  of 
constraint  coefficients  for  “factored”  constraints  with  right-hand  side  vector  b.  and  -/  the 
negative  of  the  identity  matrix.  In  this  general  development,  we  refer  to  the  F-type  const rait.ts 
as  “factored”  only  to  distinguish  them  from  the  “explicit”  E-type  constraints,  and  assume  notli- 
ing  al>out  their  structure.  !Vot  until  our  specializations  later  will  we  impose  special  structure  on 
F,  and  the  structures  we  will  consider  may  permit  the  representation  of  the  7-'-'.ype  constraints 
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without,  the  inversion  of  a  matrix.  VVe  will  set>  that  thi.s  approach  is  centercHl  around  handling 
the  part  of  the  basis  corresponding  to  the  E-type  constraints  explicitly  while  factoring  the  por¬ 
tion  of  the  basis  corresponding  to  the  E-type  constraints,  rhe  notation  is  chosen  to  suggest 
this  idea. 

Recall  that  a  basis  for  the  primal  algorithm  consists  of  ji  linearly  indejjendent  rows  froisi 
the  constraint  matrix  when  it  is  assumed  to  include  both  structural  (explicit  and  factored)  and 
nonnegativity  constraints.  Assume  that  the  current  row  basis  consists  of  k  rows  from  E,  /  rows 
from  E  and  (n  -  (k  + 1))  rows  from  — /.  Repeating  our  notation 


tc .  i  «-(*  +  /) 


B 


( 

A\\ 

\  0 


\ 

A\2 


}k  +  l 

}  ti  -  {k  +  /)  , 


where  [Aj]  A12)  includes  all  basic  .structural  rows,  from  both  E  and  E. 

We  will  ultimately  be  interested  in  i.solating  the  effect  of  each  ty|)e  of  structural  constraint 
algebraically  in  the  factored  tableau,  and  thus  we  require  greater  resolution  in  our  ffutonsl 
basis.  Introducing  obvious  notation,  we  have 


( 


[All  ^12]  = 

where  the  kernel  of  dimension  {k  -t- 1)  is  given  by 


£11 

E|2 

£1.3 

'l.  £21 

£22 

£23 

j 

[^1.! 


Ell  ^\2 
E21  E22 


Because  An  is  a  basis  foi  it  follows  that  it  is  always  possible  to  identify  among  the 

columns  of  [E21  E22]  a  nonsingular  submatrix  E22  of  dimension  1.  since  otherwise  the  rank  of 
[E21  E22]  i.s  at  most  (/  -  1)  and  thus  the  rank  of  An  is  at  most  {k  +  I  —  1),  or  equivalently 
the  rank  of  B  is  at  most  (n  -  1).  We  will  later  see  that  one  of  the  important  implementation 
challenges  is  the  task  of  efficiently  managing  the  structure  and  nonsingularity  of  En- 
The  full  factored  row  basis  is  then 


B  - 


0) 

(;» 

UJj) 


( 

Eu 

E21 

\  0 


i 

Ei2 

E22 

0 


n{k^l) 

Ei3 

/^23 

-/ 


}k 

)l 

}  n  -  {k  1) 


Introducing  the  notation 
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■^11  =  -Kn  -  E]2f'22^ F'-zt 

i4)3  =  Ei3  -  E\2F22  F'zz 

where  A\\  is  the  Schur  complement,  or  Gauss  transform,  of  F22  in  A\\  (e.g.,  Golub  and  Van 
Loan  [1983]),  we  can  write  its  inverse  as 


-F2VF2l/ir]‘  (/  +  F2-2‘F2,iri‘ £^12)^22'  F22'(F23  -  F2,ir,'/il3) 
0  0  -/ 


Grouping  the  coefficients  of  the  nonbasic  constraints  and  applying  the  same  column  ordering 
yields 


D  = 


{k+i) 


(i) 

-I 

0 

0 

)k 

(ii) 

0 

-r 

0 

V 

(Hi) 

F31 

F32 

F33 

}  Tn  -  k 

(iv) 

^  F4, 

F42 

F43  ) 

}p-l 

The  principal  part  of  the  factored  tableau  is  DP,  where  P  =  -B~'  is  the  conjugate  row 
basis.  With  the  additional  notation 


^31  =  F31  —  E22F22  F2\ 

■^33  =  F33  —  F32F22'F23 

Al  =  F41  -  F42F22'F21 

F43  =  F43  -  F42F22' F2:i 

the  principal  part  of  the  factored  tableau  is 


0) 


0» 


U}}) 


(<) 

f 

—  y4j|’  E12F jj2* 

^n'  ■^1-3 

DP  =  <"> 

(/  +  F22'F21^n'j5l2)F22' 

F22HF23  -  F2iy4,  j'/4i3) 

(Hi) 

-/iaiyiif,’ 

(A3\A^^  E]2  -  Ez2)F22 

^.33  -  A31>4,‘,'  /ii3 

(iv) 

1  -Ai/ir' 

(f^A^lEu  ~  /'V2)/-22‘ 

Fa^  -  Fii/i)|'^i3  y 
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Partitioning  w  =  (w\,W2,W3,Wi),  —  (rJ^,rT)  and  6^  =  (bj.b^)  and  introducing  the 

notation 


W-) 

W3 

62 

r\ 

f-z 


the  complete  factored  tableau  is 


wz  -  W]F22  Fz\ 
U)3  -  1^2  F22'  F23 
62  -  F41  F22'f»I 
ri  —  E\\F{2b\ 


=  r‘2 


-  F2iF22  1>1 


-^11 

-Ar/£i2F2o' 

A7/ Ai3 

-F22'F2_,Ar/ 

iI  +  Fj2'F2iAilEu)Fi2' 

F-zz'iFz'i  -  F2|/ij'|'A 

(A31A1/F12  -  F32)F22* 

A33  -  43ii4,,'Ai3 

(F41^n*^12  -  ^42)F22' 

F43  —  F41  Aj/j4i3 

.  -W2A1I 

(W2A{iEi2  ~  Wi)F22' 

W3  -  u}2Aj‘|'/4|3 

;  - 1. 


^ii'  ’1 


,3)  F22’(/'1  - 

f2  -  i43|/i  |  ,'7"i 
by  -  Ai/ij'j'n 

We  see  that  with  knowledge  of  the  current  factorization,  we  can  construct  the  entire  talileau 
from  •^22* >  -^ll*  and  the  original  problem  data.  The  dimension  of  ^22'  is  equal  to  the  number  of 
F-type  constraints  that  are  currently  basic,  and  thus  can  l>e  at  most  p.  I'he  dimension  of  A^\' 
is  equal  to  the  number  of  F-type  constraints  that  are  currently  basic,  and  thus  cannot  exctHtd 
m.  We  call  the  explicit  transformation  kernel  and  F22'  the  factored  trnns formation 
kernel. 


7  Factored  Column  and  Row  Generation 

Consider  generation  of  column  r  from  the  principal  part  of  the  tableau  DP.  Rewriting  in  a 
manner  that  highlights  our  intentions 


U) 

Ui) 

0;;) 

(i) 

^  [^n'i 

(-Ar,'Fi2F22'] 

-  -  \ 

(ii) 

{-F22'f2.[  ]} 

{^T2'-^22'^2i[ 

]} 

{F22'F23-F22'F2,[  ]} 

(Hi) 

-Ezxl  1  -  F32{  } 

-F31I  1  —  F32{ 

} 

E33,  -  Fail  1  -  F32{  } 

{iv) 

^-F4,[  ]-F42{  } 

]-F42{ 

1 

F43  -  F4i[  ]  -  £42!  }  j 

where  A\%  is  defined  as  before,  and  the  brackets  “{ j”  and  “{  }”  contain  terms  common  to  hut 
displayed  only  once  for  each  column. 

Assume  we  want  to  place  column  c  into  work  array  r.  We  partition  ~  conformal)ly  as 
refer  similarly  to  components  of  unit  vector  and  employ  a  (m  +  p)- 
vector  work  array  z.  The  notation  ”  denotes  simple  assignment,  indicates  that  a  set  of 
factored  equations  must  be  solved,  and  “=>”  explains  the  corresponding  result. 

If  column  r  is  in  (j). 


n 


and  then 


-  < - F2\[zi\ 

F22Z2  =  z  Z2  ^  {-F.^2  F2i{z]]}, 


then 


Z3  * - £3i[-=i)  -  £'32(22}, 


and  finally 


24  < - £41(21!  -  £42(22}; 

if  c  is  in  (jj),  solve 

£2222  = -(62)'"  =>  22- - (£22')*' 

z  *- EpZ2  2- - {E12F22Y 

2l  *—  j4|]*2  ^  2j  «—  [  — y4ii'£]2£22*P, 

then 

2  (62)"  -  £21(21] 

£2222=2  =»  22  ^  ((£22'^  -  £22'£2l(2l]}, 

and  finally 

23  < - £31(21]  -  £32(22} 

24  < - £41(21]  —  £42(22}; 

if  r  is  in  (jjj), 

2  (£23)*' 

£2222  =  2  =>■  22  ^  £22’ (£23)'^ 

2  *- (fiis)*"  -  £1222  =>  z  *-  {A]3Y 

Z]‘—A^^z  =>  21  *- [^7,'J4l3]^ 

then 


2  (£23)'^  -  £2121 

F2222  =  2  =►  22  -  (£22‘(£23r  -  £22'£2i[~'i1}, 

and  finally 


-3  {E33Y  ~  /^3l(^l)  -  -feVif-'i} 

^4  -  {F43r  -  F4i{c,)  E42{r2}. 

Similarly,  row  r  can  be  generated.  The  principal  part  of  the  tai>leau  i.s  now  viewed  as 


{-(  ]Fi2F2-,'} 


[  IF, 3  +  {  )F23 


[)P  ^  (ii)  [-F22'F2,ir.'i  {F22'  -  (  )F.2F22‘}  I  JF,3  +  {  }F23 

(Hi)  [-A3,iril  {(-F32-[  ]F,2)F,2'}  F.%3  +  [  ]F,3  +  {  }F23  , 

{iv)  [-Am,-/!  {(-F42  -  [  ]En)F2-/]  F43  +  [  ]F,3  +  {  }F23  y 

where  ^43,  and  Al  are  defined  as  before,  and  the  brackets  “[  ]”  and  “{  }”  contain  terms  common 
to  but  displayed  only  once  for  each  row. 

Assume  we  want  to  place  row  r  into  work  array  z.  We  partition  z  conformably  as  i  = 
(55,76.27),  refer  similarly  to  components  of  unit  vector  Cr,  and  employ  a  (rn  +  p)-vec  tor  work 
array  z. 


If  row  r  is  in  (t), 


2r> 


26F22  =  2 


il2  =>  -  ' - [.A,)'jrFi2 

^6  *—  {-[Af/)rFl2F22‘ }, 


and  finally 


if  r  is  in  (ii),  solve 


27  [25)^13  +  {z6}F23; 


26F22  =  -(e6)r 

Z  Z6F21 

25  *—  7 A,,' 


=>  26  « - (F^a'),. 

^  2*--(f22')rF2,_ 

25  *—  [  -  F22 '  F2 1 A  1 ,. , 


2  <—  (e6)r  -  [25lFi2 

26F22  =  2  =»  ^6  —  {(F22')r  -  [2r,lE,2F22' }, 


and  finally 


2?  ♦—  +  {i6}^23; 

if  r  is  in  (Hi), 

5--(E32)r 
^6^22  =  5 

Z  *—  (E3l)r  +  ZsF21 

rj  «—  (-^11*) 

then 

-  ^ - (£^32)r  -  [-Z5]£'i2 

56^22  =  5  fe  -  {(-(£32)r  -  [ir>jE,2)F,2' } 

and  finally 

-7  (£'33)r  +  [^5j^'13  +  {-=6}f^23; 

if  r  is  in  (iv), 

Z^-{Fi2)r 

26^22  =  -  =>  26  " - -_{Fi2)rF22 

Z  (F41_)r  -  Z^F2\  =>  Z  *—  (Al_)r 

25  *- Mr,']  =»  ir,  [-Ai^n  ]r 

then 

z  - {F/a)T  —  [25]£i2 

hF22  =  Z  =>  26  ^  {(-(F42)r  -  [2r,]£^12)F22’}, 

and  finally 

Z7  <—  (f43)r  +  [isj^iS  +  {26}F23. 

8  The  Complete  Algorithm 

The  complete  algorithm  is  described  in  terms  of  abstract  functions  which  operate  on  funda¬ 
mental  data  structures. 

Tableau  management  requires  two  index  maps:  one  yielding  the  intrinsic  coordinate  in  t  he 
principal  part  of  the  tableau  for  each  original,  extrinsic  problem  row  or  column,  and  the  ot  her 
its  inverse  map.  Intrinsic  arguments  arc  shown  in  lower  case,  and  extrinsic  in  upper  case. 
Index J^xchange(indexl,index2)  updates  these  maps  for  the  exchange  of  a  pair  of  tableau 
coordinates  index  1  and  index2. 

The  tableau  regions  are  successive  partitions  of  indices 


=>  2-6  -  --(E32)rF22' 

^  2  *  (2131  )r 

25  —  [-i3iv4,,']r 
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TABLEAU 

REGION 

ENDING 

INDEX 

CONTENTS  OF  REGION 

ii) 

MEC 

basic  Columns  solving  Explicit  rows 

(ii) 

MFC 

basic  Columns  solving  Factored  rows 

(Hi) 

MER 

nonbasic  Explicit  Rows 

(iv) 

m 

nonbasic  Factored  Rows 

0) 

NER 

basic  Explicit  Rows 

(jj) 

NFR 

basic  Factored  Rows 

iin) 

m  +  n 

nonbasic  Columns. 

Increment(endingJndex)  and  Decrement(endingJndex)  are  functions  to  modify  tht>se 
ending  indices. 

Generate_Row(row)  and  Genc»ate_CoIumn(column)  place  numeric  values  of  a  talileau 
row  or  column  in  ROWCOL(),  which  is  commensurate  with  the  tableau  dimensions. 

Using  ROWCOL(),  Update_Rim(row,col)  maintains  current  numeric  values  of  t  he  right- 
hand-side  and  bottom  row  of  the  tableau  in  RIM(). 

The  explicit  transformation  kernel,  is  operated  on  by  functions  using  ROW(’OL{): 

AddJE>-KernelJlow(ROW), 

Add_E-Kernel_Column(COLUMN), 

Delete_E-Kernel_Row(ROW), 

DeIete_E-Kernel_Column(COLUMN), 

ReplaceJE-Kernel_Row(REPLACEDJlOW,REPLACIN(JJlOW).  and 
UpdateJExplicit-TransformationJKernel,  the  pivotal  update. 

Aj"]'  can  be  represented  in  any  way  that  suits  the  implementer  and  efficiently  supports  these 
functions.  Generally,  we  End  A^i'  to  be  relatively  sparse,  and  report  here  results  using  an  array 
with  an  entry-point  for  each  inverse  row  and  column  which  accesses  a  stack  of  orthogonally- 
linked  nonzero  inverse  elements.  We  have  also  implemented  a  dense  vector- processor  version, 
and  the  design  i.ssues  for  LU-based  schemes  are  given  by  Ol.son  [1989]. 

The  factored  kernel,  F22,  Ls  manipulated  with: 

Factored_Kernel_Singular(ROW, COLUMN),  a  Boolean  function, 
Find_Ii-KerneLColunin_for_Key(ROW,COLUMN),  a  column  from  region  (i), 
Find_F-KerneLColumn_to_Remove(ROW, COLUMN),  a  column  from  region  (ii).  and 
Update_Factored_Kernel(ROW, COLUMN),  the  pivotal  update. 


ir 


The  complete  aljstract  algorithm  is: 


0.  Initialize; 

1.  select  primal  or  dual  algorithm; 

2.  Select  a  primal  (col)  or  dual  (row)  violation; 

STOP  if  current  solution  is  terminal,  or 
Generate-.Column(col)  or  Generate_Row(row); 

3.  By  ratio  test  with  R1M()  and  ROWCOL(), 
select  a  (row)  or  (col)  pivot  coordinate; 

STOP  if  current  solution  is  terminal,  or 
Generate_Row(row)  or  Generate_Column(coI); 

4.  Update  JRim(row,col), 

primary  Index-Exchange(row,col), 
perform  secondary  and  tertiary  exchanges  (Table  8-1) 
Update-Factored  _Kernel(ROW, COL) 
Update-Explicit-Transformation  -Kernel , 

Co  to  step  1 . 
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shown  above  may  be  required  to  preserve  the  factor 
and  Update-Explicit. Trauflformation-Kerael. 


Functions  for  maintaining  the  factored  kernel  vary  with  the  factorization,  and  a  good  im¬ 
plementation  will  exploit  these  differences.  How'ever,  a  general  specification  will  suffice  for  all 
factorizations  here. 

We  are  exclusively  interested  in  performing  two  fundamental  operations: 

1.  Solving  factored  equations  of  the  form: 

F22-2  =  ^2 

and 

~jF22=bl 

where  Z2  and  Z2  are  unknown  and  62  and  62  are  rational  (not  necessarily  integer),  and 

2.  Restoring  F22  to  the  desired  form  of  F22  which  makes  the  factored  equatioas  above  easy 
to  solve,  where  F22  results  from  inflicting  F22  with  a  column  exchange,  a  column  and  low 
deletion,  a  column  and  row  addition,  or  a  row  exchange. 

Factored_Kernel_Singular(LABELl,LABEL2)  predicts  whether  F22  will  be  singular  if: 

1.  LABELl  gives  a  basic  column  in  a  basic  factored  ROW  for  which  COL=LABEL2  would 
be  exchanged; 

2.  LABELl  is  a  basic  factored  ROW  solved  by  basic  column  COL=LABFjL2,  both  of  which 
would  be  removed  from  the  factored  kernel; 

3.  LABELl  is  a  nonbasic  factored  ROW,  which  would  be  added  to  the  factored  kernel  with 
column  COL=LABEL2,  or 

4.  LABELl  is  a  basic  factored  ROW  which  would  be  replaced  by  some  other  row  LABEL2. 

The  first  case,  a  column  exchange,  is  equivalent  to  asking  whether  solving  F22~2  =  62  with  i>2 
equal  to  region  (ii)  of  the  proposed  entering  column  COL,  yields  a  nonzero  term  associated  with 
the  basic  factored  row  in  Z2(row).  This  is  easy  to  answer  for  the  factorizations  we  discuss.  For 
instance.  Brown  and  McBride  [1984]  show  that  back-solving  a  “nearly-triangulated  generalized 
network”  basis  F22  only  as  far  as  ROW  will  suffice,  and  that  this  can  be  implemented  as  a 
traversal  of  the  “backpaths”  of  the  (zero,  one,  or  two)  coefficients  in  62  for  the  column  now’ 
basic  in  ROW.  If  numerical  precision  is  an  issue,  Z2(row)  can  be  computed  during  this  search 
and  tested  for  significance. 

The  second  case,  a  row  and  column  deletion,  is  trivial  becau.se  the  resulting  F22  will  always 
be  nonsingular. 

The  third  case,  a  row  and  column  addition,  can  be  answered  by  adding  ROW  and  COL  to 
F22  without  disturbing  its  desirable  special  structure,  thus  creating  an  instance  of  case  one.  For 
instance,  placing  ROW  first,  and  COL  last  in  F22  suffices  for  all  factorizations  di.scussed  here. 

The  fourth  case,  a  row  exchange,  can  be  amswered  by  applying  the  second  case,  deleting 
ROW  and  its  basic  column,  and  then  (perhaps  repeatedly)  applying  the  third  case,  adding  the 
row  with  index  COL  and  any  column  which  would  yield  a  nonsingular  F22- 

Find JF-KerneLColumn_to JFlemove(ROW,COL)  given  basic  factored  row  ROW,  rt*- 
turns  its  basic,  or  “key”  column  COL. 

Fmd_E-Kernel_Colunin_for_Key(LABEL,COL),  given  either  a  nonbasic  factored  ROW— 
LABEL,  or  a  basic  column  LABEL  solving  a  basic  factored  ROW,  searches  for  an  acceptable  ba¬ 
sic  column  COL  in  Explicit  Rows  (region  (i))  msing  Factored_KerneLSinguIar(HOW,COL). 

Update_Factored_KerneI(LABELl  ,LABEL2)  restores  F22  to  the  desired  form  of  F22. 
with  a  possible  increase  or  decrease  in  dimension.  Following  the  case-by-case  scheme  of  Fac- 
tored_Kernel_Singular,  we  pre-  and  post-process  the  factored  basis  repre,sentation  to  permit 
u.se  of  a  single,  static  factorization  update  of  conventional  design.  Olson  [1989]  pursues  this  in 
considerable  detail. 

Table  8-2  displays  supporting  data  structures  for  these  functions. 
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TABLE  8-2  Factorization  Algorithm  Data  Structures 


FACTORIZATIONS 

DATA 

STRUCTURE 

SIZE 

USE 

GUB,PN,GN 

RIM() 

m  +  n 

current  tableau  right-hand  .side. 

ROWCOLO 

m  +n 

bottom  row 

current  tableau  row  and  column 

MSKRCO 

m  +  n 

logical  mask  true  for  corresponding 

LQRC() 

m  +  n 

nonzero  in  ROWCOL(),  false  otherwise 
LIFO  queues  of  nonzero  row  and 

KEY() 

m  -1-  n 

column  coordinates  in  ROWC'OL() 
basic  column  in  basic  factored  row, 

PN,GN 

WORK() 

P 

and  vice  versa 

values  for  62  ^2 

MSKWKO 

P 

in  basic  factored  equations 
logical  mask  for  nonzero  row  and 

LQWK{) 

P 

column  coordinates  in  WORK() 

LIFO  queue  of  nonzero  coordinates 

PO() 

P 

in  WORKO 

next  basic  factored  row  in  pre-order 

P() 

P 

off-diagonal  row  with 

D() 

P 

nonzero  factored  coefficient 
depth,  remaining  back-substitution 

GN 

VGN() 

P 

path  length  in  factored  component 

generalized  network  cycle  factors 

JMUL() 

n 

ratio  of  generalized  network 

coefficients  in  each  column 


Generalized  upper  bound  (GUB),  pure  network  (PN),  and  generalized  network 
(GN)  factorizations  respectively  require  more  data  structures  to  support  kernel 
factorization.  For  each  factorization,  F22  is  maintained  in  some  partial  ordering 
of  rows,  and  of  columns  —  a  signed  identity  for  GUB,  upper-triangular  for  PN 
(e.g.,  Bradley.Brown,  and  Graves  [1977]),  and  nearly-upper-triangular  for  GN  (e.g,, 
Brown  and  McBride  (19841).  Direct  solution  of  factored  kernel  equations  Fnz,  =  62 
and  zj Fi2  =  bj  performed  alternately  using  the  data  structures  shown. 
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9  Computational  Experience 

The  factorization  methods  introduced  here  have  been  implemented  and  used  to  solve  a  variety 
of  existing  models  provided  by  our  colleagues.  With  their  help,  we  have  extracted  suggested 
problem  instances  from  a  diversity  of  decision  support  systems  on  host  computers  ranging 
from  mainframes  to  micro-computers.  Because  our  goal  is  to  test  factorization  technologv  in 
isolation,  the  results  reported  here  are  achieved  without  benefit  of  any  model-specific  knowhxjge 
or  tuning. 

However,  we  also  seek  to  develop  effective  modeling  tools  for  customized  use  in  developing 
and  refining  new  models.  To  this  end,  we  have  greatly  benefited  from  the  experience  and  advit  e 
of  our  colleagues,  and  we  include  some  discussion  of  modeler  guidance  and  insight  along  with 
the  numerical  results. 

Each  model  Ls  introduced  below  by  a  short  .synopsis.  Multiple  instances  of  some  models 
are  reported  where  diversity  of  size,  structure,  and  taxonomy  have  proven  interesting  to  the 
modelers.  For  those  models  that  employ  nonlinear,  mixed  integer,  or  decomposition  features, 
we  report  solution  stati.stics  for  the  initial  linear  program. 

•  GTE  The  seven  Telephone  Operating  Ctompanies  within  GTE  have  adopted  an  integrated  business  system 
called  Capital  Program  Management  System  (CPMS)  to  guide  their  3  billion  dollar  per  year  capital  planning. 
The  system  includes  a  large  scale  mixed  integer  programming  optimization  system  that  optimizes  the  critical 
economic  tradeoffs  between  maximizing  the  long-term  budget  value  of  the  firm’s  equity  and  satisfying  shorter- 
term  financial  constraints,  resource  limitations  and  service  objectives.  Investment  opportunities  for  the  next 
5  years  are  modeled  as  0-1  variables  with  alternative  implementations  for  each.  The  objective  is  to  maximize 
the  net  present  value  of  the  capital  portfolio.  There  are  financial  constraunts  on  capital,  internally  generated 
funds,  net  income  to  common,  and  limits  on  resources  such  as  labor  hours,  lines  installed,  etc.  There  are  also 
constraints  that  enforce  logical  relationships  among  opportunities  (such  as,  if  choose  A  theii  must  choose  B). 
See  Bradley  (1986|. 

•  INVEST  Capital  allocation  and  project  selection  for  Mobil  Oil  Corporation  are  modeled  as  a  two-stage 
multi-yeeir  nonlinear  capital  budgeting  problem  with  over  40,000  integer  variables.  A  master  problem  allt>- 
cates  capital  among  markets  over  a  multi-year  horizon  considering  the  estimated  nonlinear  effects  on  .sales 
of  concentrated  marketing  investments.  The  instance  reported  here  is  a  mixed-integer  linear  program  sub¬ 
problem  of  the  two-stage  model  which,  given  these  annual  cap'tal  expenditure  limits  for  a  market,  selects 
peu-ticuleir  alternate  investments.  Such  subproblems  are  easy  to  solve,  and  optimality  is  achieved  with  a  single 
iteration  of  the  nonlinear  master  problem.  See  Harrison,  Bradley  and  Brown  [1989]. 

•  TANKER  A  crude  oil  tanker  scheduling  problem  faced  by  a  major  oil  company  is  solved  using  an  elastic 
set  partitioning  model.  The  model  takes  into  account  all  fleet  cost  components,  including  the  opportunity 
cost  of  ship  time,  port  and  canal  charges,  demurrage  and  bunker  fuel.  The  model  determines  optimal  speeds 
for  the  ships  and  the  best  routing  of  ballast  (empty)  legs,  as  well  as  which  cargoes  to  load  on  controlled  ships 
and  which  to  spot  charter.  All  feasible  schedules  are  generated  the  cost  of  each  is  accurately  determined 
and  the  best  set  of  schedules  is  selected.  See  Brown,  Graves  and  Ronen  [1987). 

•  HFDF  A  large-sctJe  elastic  set  partitioning  model  used  to  assign  frequencies  for  a  network  of  high  frequency 
direction  finding  receivers.  See  Brown,  Drake,  Marsh,  and  Washburn  [1990|, 

•  CAS  A  multi-time  period  strategic  model  for  use  by  natural  gas  utilities  for  determining  optimal  contract 
levels  for  gas  purchase,  storage  and  transmission.  An  underlying  generalized  network  flow  model  represents 
gas  being  bought,  stored,  shipped  and  consumed  over  a  multi-year  time  horizon,  typically  at  a  monthly  level 
of  detail.  Constraints  and  variables  are  added  to  handle  variable  maximum  and  minimum  purcheise  levels, 
variable  leased  or  constructed  storage  and  variable  transmission  capacities.  An  integrated  parallel  model 
incorporates  the  peak  requirements  necessary  on  some  days  during  cold  winter  months.  This  model  has  been 
used  by  a  number  of  utilities  including  Southwest  Gas  C^orporation  and  Questar  Pipeline  (Corporation  to  plan 
operations  and  to  justify  such  plans  to  regulatory  agencies.  See  Avery,  Brown,  Rosenkranz  and  Wood  [1992] 

•  KELLOGG  A  multi-time  period,  multi-plant  production/invcntory/tran.sshipnient  linear  program  for  Kellogg 
cereals.  The  model  guides  weekly  processing,  packaging  and  shipping  decisions.  Production  consists  of  two 
stages:  processing  lines  produce  basic  products  which  are  then  packaged  on  packaging  lines  into  different- 
sized  containers  to  yield  finished  products.  Processing  lines  produce  a  subset  of  the  basic  Troducts  and  have 
limited  capacity  with  overtime  charges  for  weekend  shifts.  Packaging  lines  are  analogous.  In-house  inventory 
rapacity  is  limited  although  outside  storage  is  available  at  additional  cost.  Inter-plant  shipments  of  finished 
products  are  made  by  rail  or  truck.  See  Wo<k1  jl989]. 
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•  ODS  A  commonly  occurring  problem  in  distribution  system  design  is  the  optima!  location  o!  intermediate 
distribution  facilities  between  plants  and  customers.  A  nuilticoiiimodity  capacitated  single- pcriiKl  version 
of  this  problem  is  formulated  as  a  mixed  integer  linear  program  A  stjiulion  technique  bas<;d  on  Menders 
Decomposition  is  developed.  . . ,  An  essentially  optimal  solution  is  found  and  proven  with  a  surprisingiv  small 
number  of  Benders  cuts.  See  Geoffrion  and  Graves  [1974].  The  instances  reportwl  here  are  deroniposition 
master  problems. 

•  TAM  The  annual  decision  on  how  much  the  Air  Force  should  spend  on  aircraft  and  on  munitions  is  of  great 
interest  to  many  people.  How  the  Air  Force  stalF  develops  information  to  support  the  decision  has  changed 
over  the  years.  Currently,  a  linear  program  is  being  used  by  the  Air  Force  Center  for  Studies  and  Analysis 
and  is  being  tested  by  the  Munitions  Division  of  the  Plans  and  Operations  Directorate  (AF/XOXFM)  for 
munitions  tradeoff  analysis.  The  LP  uses  existing  data  and  estimates  of  (1)  aircraft  and  munition  effectiveness, 
(2)  target  value,  (3)  attrition,  (4)  aircraft  and  munition  costs,  and  (5)  existing  inventories  of  aircraft  and 
munitions.  Other  factors  considered  are  weather  and  length  of  the  conflict.  See  Might  [1987]  and  Jackson 
(1989). 

•  PHOENIX  A  planning  model  for  the  multi-year,  multi-billion  dollar  modernization  of  the  F.  S.  Army's  aging 
helicopter  fleet.  The  mixed  integer  linear  program  employs  a  multi-product  production/iiiventory  formulation 
with  aged  inventory.  Goal  constraints  attempt  to  enforce  fleet  size,  maximum  age,  and  teclinolivgy  goals  for 
each  year  and  each  of  four  aircraft  missions,  while  also  keeping  expenditures  within  upper  and  lower  limits 
Additionally,  combinatorial  constraints  and  variables  handle  production  line  startup  and  shutdown  costs, 
minimum  and  maximum  production  levels  and  requirements  linking  certain  production  lines  See  ('Icmence. 
Teufert,  Brown  and  Wood  [1988]  and  Brow.n,  Clemence,  Teufert,  and  Wood  [I990j. 

•  EA6B  Configures  jammers  of  hostile  radar  on  an  EA-6B  “Prowler”  .Xaval  electronic  warfare  airrrafi.  See 
Sterling  [1990). 

•  DEC  Digital  Equipment  Corporation  uses  this  model  to  determine  worldwide  manufacturing  and  distribiiiion 
strategy  for  new  products.  This  mixed  integer,  linear  program  suggests  a  production,  distribution,  and  vendor 
network  which  minimizes  cost  and/or  cumulative  cycle  times  subject  to  constraints  on  estimated  demand, 
local  content,  and  joint  capacity,  over  multiple  products,  echelons,  and  time  periods.  Ck>st  fact^'rs  include 
fixed  and  variable  production  charges,  distribution  via  multiple  modes,  taxes,  duties  and  duty  drawback,  and 
inventory  charges.  See  Harrison,  Arntzen,  and  Brown  (1992). 

•  AMMO  4H  A  four-commodity  transshipment  model  for  delivery  over  time  of  military  products  from  pro 
duction  and  storage  locations  to  overseas  locations  to  support  theater  operations  is  developed  The  model 
covers  five  physical  echelons,  including  production  plants,  .storage  depots,  ports  of  embarkation,  ports  of 
debarkation  and  geographic  field  locations.  Road,  rail,  sea  and  air  transportation  are  modeled,  and  product 
demands  are  time-phased.  Capacitation  occurs  primarily  on  sea  and  air  links,  and  as  throughput  capacities 
on  transfer  points,  requiring  replication  of  some  echelons.  The  ob  lective  of  the  model  is  to  minimize  deviation 
from  on-time  deliveries.  See  Staniec  [1984]. 

•  BUSCH  A  model  of  brewery-to-wholesaler  movements  of  beer  for  Anheuser-Busch.  The  model  also  includes 
some  packaging  decisions  and  is  essentially  a  multicommodity  flow  model  with  joint  capacity  constraints 
arising  from  loading  dock  and  inventory  capacities  as  well  as  some  managerial  requirements.  See  Brown, 
Mamer,  McBride  and  Wood  [1992].  The  instance  reported  here  is  a  small  pilot  model  for  the  full-scale 
system  with  millions  of  variables  which  is  solved  directly,  or  by  decomposition. 

•  BAB-  A  linear,  mixed-integer  multi-period  production-inventory  master  planning  model.  See  Harrison  |1992j. 

Four  implementation.s  are  compared:  “XS”  is  an  unadorned  version  of  the  X-Systein,  an 
implementation  of  the  Graves  mutual  primal-dual  method  with  its  GUB  factorization  disabled, 
while  “XS{GUB)”,  “XS(PN)”,  and  “X.S(GN)”  each  employ  the  respective  factorizations  dis- 
cus.sed  here.  To  estlablish  a  frame  of  reference,  performance  of  these  implementations  is  com¬ 
pared  with  two  well-known  commercial  solvens:  IBM’s  Optimization  Subroutine  Library  “OSL” 
(Relea.se  2  [1991]),  and  two  versions  of  “CPLEX”  (Version  1.2  [1990]  and  Version  2.0  [1992]). 

Ideally,  one  would  develop  four  equivalent  formulations  of  each  model,  each  customizfHl  for 
its  particular  solver  with  the  goal  of  inducing  a  large  factored  row  set  of  the  appropriate  type. 
This  approach  is  a  consistent  theme  in  the  literature  dealing  with  specialized  algorithms  and 
one  that  we  strongly  endorse.  Alternate  formulations  of  a  model  are  often  available,  and  it 
.seems  scn-sible  to  choose  one  that  exploits  as  much  as  possible  the  strengths  of  the  solver. 


However,  all  of  the  models  used  here  are  “off-the-shelf”  in  the  sense  t  hat  the\  were  deveiop<Nl 
at  various  times  by  various  modelers,  and  alternate  formulations  are  imj)rai  tiral)le.  Thus,  the 
approach  is  to  preserve  a  single,  unfactored  representation  of  each  model,  and  attempt  to 
identify  favorable  row  structures  through  the  use  of  heuristics.  I'he  proctsiure  is  basts!  on 
the  work  of  Brown  and  Thomen  [1980),  Brown  and  Wright  [198:jj  and  Brown,  .McBride  ami 
Wood  (1985].  The  heuristics  are  greedy  and  myopic  in  the  sense  that  tfiev  initially  consider  tlu' 
entire  row  set  of  the  problem,  and  discard  one  row  at  a  time  without  backtracking  until  the 
remaining  set  satisfies  the  desired  row  factorization.  This  can  be  expected  to  confound  or  dt^'roy 
structure  introduced  by  the  modeler.  Although  the  automatic  factorization  implementation  htis 
options  to  accept  modeler  guidance,  the  methods  are  compared  here  without  this  subjf'ctive 
complication.  While  these  model-naive  experiments  yield  interesting  and  uw'ful  oltservat ions 
about  the  implementations,  they  suffer  for  lack  of  guidance  by  a  skilled  modeler. 
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Table  9-1  shows  the  important  structural  information  concerning  the  model  instam(?s  to  l>e 
solved. 


TABLE  9-1  Problem  Dimensions 


n 

rn 

Vai’B/% 

Pfn  1% 

NZEL 

CTE 

6,624 

960 

909/95 

909/95 

922/96 

58 

INVEST 

11,989 

1 ,338 

941/70 

1,101/82 

1,168/87 

.33 

TANKER 

7,598 

83 

32/39 

32/39 

66/80 

31 

HFDF 

10,548 

61 

31/51 

31/51 

:\2jTi2 

]89 

GAS  PN  A 

27,884 

6,848 

4,345/63 

5,9:14/87 

5.976/87 

37 

GAS  PN  C 

15,362 

3,794 

2,658/70 

:i,418/90 

:Ld20/90 

20 

GAS  PN  E 

5,102 

1,184 

434/37 

877/74 

88.3/75 

7 

GAS  GN  A 

27,884 

6,848 

4,484/65 

5,142/75 

5.976/87 

37 

GAS  GN  C 

15,362 

3,794 

2,664/70 

:J, 084/81 

3,420/90 

20 

KELLOGG  2 

17,841 

3,818 

1,265/33 

2,578/68 

2..596/68 

:i5 

KELLOGG  3 

27,490 

5,727 

2,295/40 

3,867/68 

3,892/68 

54 

KELLOGG  4 

37,139 

7,636 

2,428/32 

5,156/68 

5,188/68 

74 

KELLOGG  5 

46,788 

9,545 

3,388/36 

6,445/68 

6,484/68 

93 

ODS  1 

11,568 

3,023 

528/17 

540/18 

558/18 

21 

ODS  3 

23,993 

594 

490/82 

490/82 

490/82 

68 

TAM  5 

10,531 

438 

102/23 

1.32/;i0 

162/37 

94 

TAM  8 

6,104 

420 

118/28 

154/37 

196/47 

49 

TAM  12 

17,793 

629 

177/28 

231 /37 

294/47 

165 

PHOENIX  10 

6,884 

1,618 

206/13 

220/14 

1,15.3/71 

14 

PHOENIX  30 

17,212 

4,305 

293/07 

303/07 

:{,604/84 

48 

EA6B 

12,247 

2,978 

1,610/54 

2,921/98 

2,921/98 

16 

DEC 

14,518 

2,171 

677/31 

677/31 

1 ,088/50 

24 

AMMO  4H 

83,497 

13,963 

6,874/49 

12,892/92 

12,892/92 

129 

BUSCH  4 

7,997 

1,248 

649/52 

1,140/91 

1,148/92 

15 

BAR 

49,032 

7,446 

2,712/36 

4,575/61 

5,1.34/69 

102 

For  each  problem,  the  total  number  of  structural  variables  is  n,  structural  con¬ 
straints  771,  GUB  rows  found  by  the  identification  heuristic  pccc,  pure  network 
rows  Pf'^,  generalized  network  rows  pcN,  and  thousands  of  nonzero  technological 
coefficients  NZEL.  For  example,  GTE  has  58  thousand  nonzero  technological  co¬ 
efficients  for  6,624  variables;  GTE  can  be  viewed  as  having  960  explicit  and  no 
factored  rows,  or  as  909/95%  GUB-factored  and  960  —  909  =  51  explicit  rows,  as 
909  PN-factored  rows,  or  as  922  GN-factore<l  and  38  explicit  rows. 


Table  9-2  displays  solution  times  for  the  sample  problems.  "I'hese  C'Pl'  times  exclude  initial 
problem  input,  factored  row  identification,  and  final  output  -  on  average  about  0.2  second  per 
problem. 


TABLE  9-2  Solution  Seconds 


AMDAHL  ri99r,-700 


486/:f3MHz  PC 


none 

X-Systeni 
GUB  PN 

CN 

IBM 

OSL 

XS 

CN 

CPLEX 

1.2  2.0 

CTE 

9 

8 

8 

8 

43 

41 

65 

58 

INVEST 

4 

5 

4 

4 

23 

24 

52 

49 

TANKER 

8 

10 

10 

8 

6 

48 

8 

9 

HFDF 

100 

101 

101 

111 

40 

462 

581 

542 

CAS  PN  A 

2,376 

778 

50 

57 

291 

321 

1,714 

1,221 

CAS  PN  C 

4 

4 

3 

3 

79 

26 

185 

245 

CAS  PN  E 

72 

33 

4 

6 

13 

33 

165 

50 

CAS  CN  A 

1,115 

542 

294 

72 

312 

385 

3,532 

2.262 

GAS  GN  C 

4 

4 

3 

3 

71 

26 

186 

236 

KELLOGG  2 

3 

2 

2 

2 

69 

5 

219 

194 

KELLOGG  3 

5 

5 

5 

5 

144 

9 

513 

440 

KELLOGG  4 

48 

36 

26 

24 

500 

115 

1,274 

1,111 

KELLOGG  5 

1,122 

1,459 

320 

248 

1,210 

926 

2,615 

2,243 

ODS  1 

6 

3 

2 

5 

125 

22 

1 ,042 

414 

ODS  3 

6 

6 

6 

6 

23 

38 

58 

56 

TAM  5 

55 

60 

45 

40 

44 

231 

151 

86 

TAM  8 

24 

14 

14 

17 

21 

69 

77 

71 

TAM  12 

108 

112 

88 

106 

101 

312 

416 

378 

PHOENIX  10 

4 

2 

2 

2 

20 

10 

84 

73 

PHOENIX  30 

41 

25 

26 

9 

335 

69 

1,171 

1,2.39 

EA6B 

75 

33 

9 

9 

45 

44 

177 

244 

DEC 

189 

32 

42 

80 

71 

93 

317 

293 

AMMO  4H 

42 

41 

35 

46 

1,000 

277 

1,762 

1,597 

BUSCH  4 

5 

5 

4 

4 

26 

34 

55 

46 

BAR 

290 

212 

106 

114 

382 

480 

1 ,366 

1,222 

An  AMDAHL  5995-700  running;  under  IBM  VM/CMS/XA  with  IBM  VS  FOR 
TRAN  2.3.0  is  used  to  render  performance  in  CPlI-seconds  accurate  to  the  precision 
shown  for  the  basic  “XS”-system  using  no  factorization,  compared  with  dynamic 
factorizations  of  '‘XS(GUB)”,  pure  network  “XS(PN)",  and  generalized  network 
'‘XS(GN)”  rows.  “IBM-OSL”  shows  primal  simplex  performance  on  the  same  com¬ 
puter  of  the  IBM  Optimization  Subroutine  Library,  Release  2  [1991].  “486/.3.3” 
shows  the  clock-time  performance  of  “XS(GN)”  on  a  microcomputer  (33  MHz  Intel 
486  with  32  MB  RAM)  followed  by  that  of  "C’PLEX”  (Version  1.2)  [1990]  and  of 
“GPLEX”  (Version  2.0)  [1992]  on  the  same  microcomputer. 


The  original  formulations  of  most  of  the  test  problems  are  strongly  influenced  by  the  mod¬ 
eler’s  solution  strategy.  For  instance,  TANKER  is  endowed  with  a  CUB  structure  which  places 
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every  binary  variable  in  an  associated  set  from  which  only  one  member  can  be  chosen:  this 
maximal  GIJB  set  is  also  sought  for  its  tendency  tt)  yield  nearly-integer  solutions  to  the  linear 
program.  AMMO  4H  is  a  multicommodity  capacitateri  transshipment  problem  and  thus  is  Itest 
suited  to  a  pure  network  factorization;  it  was  originally  solved  ity  dual  dccomiiosition  rendering 
pure  network  sulr-problents.  PHOENIX  is  a  multicotnmodity  erjuipment  replacentent  mode! 
closely  following  the  generalized  network  factorization  paradigm. 

The  row  structures  in  Table  9-2  have  been  found  without  modeler  help  by  our  automatic 
identification  heuristic,  yet  they  corroborate  the  modelers’  intentions.  TANKER  reveals  the 
same  pure  network  rows  as  the  Ol'B  rows,  or  a  generalized  network  constructwl  by  identifying 
one  additional  row  to  be  paired  with  each  Gl.'B  row.  f’HOENlX  exhibits  a  dominant  structure 
that  is  clearly  a  generalized  network. 

One  would  anticipate  the  factorization  exploiting  the  dominant  row  structure  to  win  com¬ 
putation  tests.  This  is  wrong  more  often  than  right.  Table  9-2  shows  that  the  more  genera! 
factorization  dominates  the  less  general,  wuth  few'  exceptions:  GTE’s  relatively  large  (U'B  set, 
and  the  large  pure  networks  in  GAS  PN  A,  GAS  F*N  E,  and  AMMO  4H  seem  to  satisfy  our 
prior  bias  toward  model-dictated  factorizations. 

Table  9-2  also  suggests  that  our  myopic  use  of  heuristics  to  automatically  identify  factorfHl 
structure  has  its  pitfalls.  In  a  number  of  problem  instances,  we  identify  significantly  larger  fac¬ 
tored  sets  with  the  more  general  factorizations,  yet  we  enjoy  little  improvement  in  computation 
times  (e.g.,  INVEST,  ODS  3,  TAM  8).  This  suggests  that  the  “quality”  of  a  row  factorization 
is  not  completely  specified  by  its  size. 

We  have  pursued  this  notion  by  inviting  some  of  the  modelers  to  guide  our  identificatiim 
heuristic  to  precisely  the  row  .sets  they  intended.  Some  of  the  results  have  been  striking:  Wood 
[1989]  reports  significant  improvements  for  problems  in  the  GAS  system. 

A  few'  of  our  corresponding  modelers  have  had  the  opportunity  to  build  models  from  scratch 
with  a  particular  factorization  in  mind.  This  admits  model  coercion  and  a  w'ide  range  of  well- 
know'n  reformulation  methods  w'hich  we  think  can  materially  change  both  the  size  and  quality 
of  the  result.  Their  early  reports  show  promi.se.  Among  the  models  discu.ssed  here,  the  GAS 
and  KELLOGG  .systems  have  been  sul)sequently  reformulated  to  enhance  generalized  networks. 
GTE  has  been  re-engineered  to  further  accentuate  its  dominant  Gl'B  set,  and  TANKER-like 
and  many  ODS  models  have  been  moved  to  a  micro-computer;  all  these  models  are  now  larger, 
but  much  easier  to  solve. 

It  is  surprising  and  encouraging  that  the  transition  to  more  general  factorizations  .seldom 
degrades  performance  much,  even  when  few  additional  factored  rows  are  won  in'  the  increase<l 
generality.  This  contradicts  popular  folklore  that  the  more  general  factorizat  ions  demand  sulv 
stantial,  if  not  overwhelming  increases  in  the  resulting  sizes  of  the  factoreri  structures.  In 
fact,  computational  testing  reported  by  others  has  usually  been  limited  to  iiiodels  in  which 
the  number  of  explicit  rows  is  in  the  range  of  one  to  twenty  (e.g.,  Ghon  and  .Saigal  [1977] . 
Glover,  Karney,  Klingman  and  Russell  [1978],  Glover  and  Klingman  [1981]).  Our  results  are 
all  the  more  remarkable  given  the  lack  of  guidance  frotn  t  he  modeler  for  tlie  “int ('tided’'  row 
factorization. 


27 


Table  9-3  shows  the  maximum  size  of  the  v4|i'  in  terms  of  its  nonzero  elements. 

TABLE  9-3  Maximum  Number  of  Elements  in  Explicit  Transformation  Kernel 


XS 

XS(CUB) 

XS(PN) 

XS(CN) 

CTE 

13 

0 

0 

0 

INVEST 

9 

1 

1 

1 

TANKER 

1 

1 

1 

0 

HFDF 

3 

1 

1 

1 

CAS  PN  A 

1,550 

891 

30 

54 

CAS  PN  C 

18 

5 

0 

0 

CAS  PN  E 

263 

no 

i 

5 

CAS  GN  A 

1,470 

758 

340 

59 

CAS  CN  C 

19 

/ 

1 

0 

KELLOCC  2 

6 

2 

0 

0 

KELLOGG  3 

9 

4 

0 

0 

KELLOCC  4 

79 

41 

10 

11 

KELLOCC  5 

1,299 

1,015 

138 

128 

ODS  1 

9 

1 

1 

5 

ODS  3 

0 

0 

0 

0 

TAM  5 

28 

22 

15 

18 

TAM  8 

20 

15 

8 

10 

TAM  12 

38 

42 

16 

21 

PHOENIX  10 

32 

21 

21 

1 

PHOENIX  30 

169 

157 

185 

1 

EA6B 

1,155 

270 

0 

0 

DEC 

218 

39 

49 

46 

AMMO  4H 

235 

92 

0 

0 

BUSCH  4 

10 

5 

0 

0 

BAR 

290 

163 

70 

41 

The  number  of  nonzero  elements  (in  nearest  thousands)  in  the  explicit  transforma¬ 
tion  kernel  gives  some  indication  of  how  much  information  is  not  captured  by 
factorization,  as  well  as  an  idea  of  relative  storage  requirements. 

We  see  that  the  maximum  size  of  the  explicit  transformation  kernel  tends  to  decrease  as  the 
generality  of  the  factorization  increases.  Recalling  the  definition  of  the  explicit  transformation 
kernel, 

A,V  =  (^11  '  . 

this  trend  is  as  we  would  expect.  Each  potentially  binding  explicit  row  which  can  be  converted 
to  a  factored  row  reduces  the  likely  size  of  A]/.  Also,  the  density  of  the  term  -EnF-u  f''2\ 
generally  increases  with  the  density  of  F'n-  For  F22  k-hy-k,  the  number  of  nonzeros  in 
for  the  CUB  factorization  is  fc,  for  pure  networks  perhaps  as  large  y,  and  for  generalized 
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networks  as  large  as  fc*.  There  are  some  exceptions  to  this  trend  in  Table  9-d.  (‘specially  in 
model  instances  in  which  the  size  of  the  (C!N)  factored  row  set  is  not  significant h  larger  than 
that  of  (PN)-  This  is  because  for  a  given  factored  kernel  /-Vj.  the  exact  (P.\)  representation  of 
F.y,^  is  generally  more  sparse  than  that  of  the  less  exact  floating-point  representation  of  ((IN). 
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It  is  usually  the  case  that  many  constraints  are  not  binding  at  optimality,  a«  can  i>e  st'en  m 
I'able  9-1. 


TABLE  9-4  Binding  Explicit  Constraints  at  Optimality 


Binding/% 

GUB/% 

PN7% 

GN7% 

CTE 

554/58 

20/02 

20/02 

18/02 

INVEST 

763/57 

201/15 

195/15 

162/12 

TANKER 

51/61 

60/72 

39/47 

9/11 

HFDF 

58/95 

29/48 

29/48 

28/46 

GAS  PN  A 

3,068/45 

1,953/29 

299/04 

349/05 

GAS  PN  C 

2,324/61 

850/22 

68/02 

90/02 

GAS  PN  E 

901 /76 

573/^8 

90/08 

79/07 

GAS  GN  A 

3,064/45 

1,8.54/27 

1,1.55/17 

3.59/05 

GAS  GN  C 

2,32.3/61 

861/2.3 

402/11 

89/02 

KELLOGG  2 

1,950/51 

1,015/27 

92/02 

118/0.3 

KELLOGG  3 

2,942/51 

1,368/24 

148/03 

18.3/03 

KELLOGG  4 

4,1.36/54 

2,491/33 

391/05 

452/06 

KELLOGG  5 

5,270/55 

2,305/24 

.592/06 

617/06 

ODS  1 

361/12 

83/03 

76/0.3 

117/04 

ODS  3 

410/69 

0/00 

0/00 

0/00 

TAM  5 

264/60 

227/52 

168/38 

186/42 

TAM  8 

279/66 

240/57 

142/.34 

175/42 

TAM  12 

412/66 

352/53 

205/.3.3 

251 /40 

PHOENIX  10 

1,098/68 

1,096/68 

1,090/67 

80/05 

PHOENIX  30 

3,298/77 

3,2.36/75 

3,2.39/75 

110/0.3 

EA6B 

2,9.39/99 

1,. 3.34/45 

26/01 

27/01 

DEC 

1,328/61 

7.36/34 

726/33 

421/19 

AMMO  4H 

6,686/46 

3,15.3/22 

13/00 

14/00 

BUSCH  4 

840/67 

481 /,39 

5/00 

8/01 

BAR 

4,687/63 

.3,2.50/44 

1,895/25 

1,4.50/19 

Not  all  constraints  are  binding  at  optimality.  The  first  column  lists  the  number 
of  binding  explicit  constraints  and  expresses  this  as  a  percentage  of  all  constraints: 
the  following  columns  display  the  number  of  binding  explicit  constraiiits  and  their 
corresponding  percentages  under  the  alternate  factorizations. 


A  distinguishing  feature  of  dynamic  factorization  is  the  ability  to  limit  attention  to  binding 
constraints,  handling  binding  factored  constraints  with  great  efficiency,  and  working  with  a 
relatively  small  number  of  Irinding  explicit  constraints.  Explicit  binding  constraints  on  the 
order  of  a  few  thousand,  or  less,  are  quite  manageable.  This  is  well  beyond  the  size  of  previously 
reported  implementations. 
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10  Conclusions 


F*reviou.s  research  by  others  generally  suggests  that  spt>cialized  algorithms  such  as  tlujse  pre¬ 
sented  here  are  useful  only  when  the  factored  structure  completely  dominates.  There  are  even 
reports  of  algorithms  for  solving  problems  having  a  single  unfactoreil  (explicit)  constraint  (Hultz 
and  Klingman  [1978],  Klingman  and  Russell  [1978]).  When  implementations  have  be«-n  r<*- 
ported,  problem  suites  have  Ijeen  limited  to  instances  having  a  very  small  number  of  explit  it 
constraints,  typically  in  the  range  from  one  to  twenty  (Chen  and  Saigal  [1977],  Clover.  Karnoy, 
Klingman  and  Russell  [1978],  Clover  and  Klingman  [1981]}.  The  cttnserisus  seems  t<j  be  that 
such  algorithms  are  quite  delicate,  and  deserve  to  be  viewed  as  specializctl  algorithms,  useful 
only  for  solving  very  special  problem  instances. 

We  refute  this  view.  Dynamic  factorization  is  competitive  with  commercial-tiuaiity  opti¬ 
mization  systems  on  every  model  instance  we  have  tested. 

The  development  here  stres.ses  the  similarities  among  the  algorithms  and  the  nat  ural  exten¬ 
sions  leading  from  one  to  the  next.  This  is  in  contrast  to  the  development  report e«l  for  similar, 
non-dynamic  algorithms  (e.g.,  Dantzig  and  Van  Slyke  []967j,  Klingman  and  Russell  [1978]. 
and  Hultz  and  Klingman  [1978])  in  which  the  specifics  of  the  individual  algorithm  olxscure  the 
generality  of  the  approach.  The  conceptual  difference  between  our  algorithms  is  sp<>n  to  be 
largely  Isolated  to  the  structure  of  a  single  algebraic  entity,  the  factored  kernel.  By  abstracting 
the  structure  of  the  factored  kernel  and  concentrating  on  the  genera!  algorithm  design,  the 
versatility  and  flexibility  of  this  approach  is  clarified. 

The  algorithmic  development  leads  directly  to  an  implementation.  The  resulting  software 
suite  exhibits  a  “.single  sy.stem  image”.  The  modularity  of  the  algorithm  allows  the  definition 
of  an  “abstract  data  type”  (see,  e.g.,  A.ho,  Hopcroft  aird  L'llman  [1974])  which  isolates  the  data 
structures  and  update  procedures  for  the  factored  kernel  from  the  rest  of  the  implementation. 
Each  factorization  is  .seamlessly  integrated  within  the  system  design. 

The  early  1980s  produced  a  great  deal  of  research  on  automatic  identification  of  special 
structure  in  LP  models  (see,  e.g.,  Clunawardane  and  Schrage  [1977],  Cllover  [1980],  Schrage 
[1981],  Brown,  McBride  and  Wood  [1985],  and  Bixby  and  Tourer  [1986]).  We  have  inrorpuratefi 
the  most  useful  of  these  ideas  into  our  implementation,  and  we  have  what  we  believe  to  be  the 
first  complete  implementation  which  supports  automatic  identification  of  factored  row  sets.  This 
capability  may  be  used  to  identify  new  factored  structure  or  to  validate  or  augment  a  modeler- 
provided  recommendation.  When  faced  with  the  choice  of  either  solving  an  unfactored  model 
instance  or  automatically  identifying  a  factored  .structure  and  then  using  the  corresponding 
.solver,  our  results  show  that  the  latter  is  nearly  always  to  be  preferred.  Modelers  have  conducted 
extensive  zidditional  computational  experimentation  with  the  X-System  not  reportwi  in  this 
paper.  These  results  suggest  that  in  addition  to  the  quantity  of  fartore<l  rows,  the  quality  of 
these  rows  influences  the  performance  of  factorization  algorithms.  W'hile  not  well  understood, 
it  is  clear  that  the  myopic  approach  of  our  heuri.stics  is  no  substitute  for  the  modelers  guidance 
in  identifying  factored  structure. 

Processing  networks  (Koene  [1982])  are  network  models  which  allow  proportional  flow  re¬ 
strictions  on  the  arcs  entering  or  leaving  .some  nodes.  One  formulation  of  such  a  model  results 
in  a  pure  or  generalized  network  .structure  with  a  set  of  complicating  columns,  t'hen  and  En- 
gquist  [1986]  propose  a  primal  partitioning  algorithm  for  solving  processing  network  problems. 
An  alternate  formulation  yields  a  pure  or  generalize<]  network  structure  with  complicating  rows: 
this  is  precisely  the  structure  dealt  with  here. 

The  multicommodity  capacitated  transshipment  prol)lem  (Mt'  l'P)  has  been  the  subject  of 
much  research  over  the  years,  and  a  numl>er  of  specialized  algorithms  have  be<'n  propo.sf>d  to 
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solve  it  (see,  e.g.,  Assad  [1978]  or  Kennington  [1978]).  The  u.sual  MCTI’  formulation  is  a  pure 
network  which  each  commodity  uses  independently  in  its  own  flt)W  model,  hut  with  side  ton- 
straints  on  the  total  common  flow  of  all  commodities  over  some  of  the  network  an  s.  The  side 
constraints  form  a  (JUB  row  set,  while  the  rest  of  MC'TP  forms  a  pure  network;  either  view 
might  be  preferred  depending  upon  size  of  the  common  network,  the  number  of  side  const  raints, 
and  the  number  of  commodities.  In  our  experience,  the  network  factorization  usually  domi¬ 
nates  the  CUB  factorization,  and  the  pure  network  factorizaticm  presented  here  is  a  powerful 
technique  for  .solving  MCTPs.  As  an  experiment,  we  customized  our  (PN)  implementation  for 
MCTP  to  exploit  the  special  structure  of  the  explicit  side  constraints.  1  his  highly-specialized 
implementation  performed  no  better  on  AMMO  4H,  and  w'e  now  believe  that  this  would  be 
true  for  most  MCTPs. 

There  are  problems  which  would  exhibit  a  large  factored  row  structure  if  not  for  a  set  of 
complicating  columns  (e.g.,  see  Brown,  McBride  and  Wood  [1987)]).  One  would  expect  the 
structure  of  the  factored  kernel  to  be  dominated  by  that  of  the  predominant  row  structure, 
with  only  occasional  complications  due  to  the  exceptional  columns.  One  might  allow  for  this 
exceptional  structure  in  the  factored  kernel  by  identifying  it  “on-the-fly”  as  the  algorithm 
progresses,  and  preserving  the  .sanctity  of  the  core  factorization.  I’hough  conceptually  simple, 
.some  iterations  of  this  algorithm  would  border  on  the  spectacular.  This  approach  may  l)e 
thought  of  as  a  hybrid  between  the  dynamic  factorization  developed  here  and  dynamif  luisis 
triangulation  methods  (see,  e.g.,  Hellerman  and  Rarick  [1971]  and  [1972].  .Saunders  [1978]  and 
McBride  [1980]). 

Dynamic  extrinsic  factorization  is  subsumed  by  the  algorithms  presente<i  in  this  paper  if 
we  activate  functions  in  the  update  analogous  to  the  secondary  exchanges  now  employed.  Es¬ 
sentially  all  that  has  to  be  done  is  ensure  that  successive  factored  components  retain  their 
stipulated  special  structure.  We  speculate  that  this  will  work  best  in  cases  where  model  struc¬ 
ture  is  amenable,  and  quite  likely  will  require  some  model-specific  customization  to  perform 
well  on  difficult  models.  We  have  limited  our  experimentation  to  those  static  extrinsic  cases 
which  we  believe  to  be  most  generally  useful. 
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